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O ; ABSTRACT 

^ . Using a cosmological N-Body simulation and a sample of re-simulatcd cluster-like 

' haloes, we study the mass loss rates of dark matter subhaloes, and interpret the mass 

function of subhaloes at redshift zero in terms of the evolution of the mass function of 
systems accreted by the main halo progenitor (hereafter called the 'unevolved subhalo 
mass function'). When expressed in terms of the ratio between the mass of the subhalo 
at the time of accretion, m^, and the present day host mass, Mq, the unevolved subhalo 
' mass function is found to be universal, in that it is independent of the mass of the host 

*~^[ halo. However, the subhalo mass function at redshift zero (hereafter called the 'evolved 

Q , subhalo mass function') clearly depends on Mq, in that more massive host haloes host 

■ more subhaloes. In order to relate the unevolved and evolved subhalo mass functions, 
C/3 I we measure the subhalo mass loss rate as a function of host mass and redshift. We 
a . find that the average, specific mass loss rate of dark matter subhaloes depends mainly 

on redshift, with only a very weak dependence on the instantaneous ratio between 
, the mass of the subhalo, nisb, and that of the host halo at that time My. In fact, to 

^ ■ good approximation, subhalo masses 'decay' exponentially, with a decay-time that is 

fT^ ' proportional to the instantaneous dynamical time of the host halo. Combined with the 

\^ . fact that more massive haloes assemble later, these results suggest a pleasingly simple 

•/^ ' picture for the evolution and mass dependence of the evolved subhalo mass function. 

^""I I Less massive host haloes accrete their subhaloes earlier, which are thus subjected to 

^SJ . mass loss for a longer time. In addition, their subhaloes are typically accreted by 

,—1 ' denser hosts, which causes an additional boost of the mass loss rate. To test the self- 

, consistency of this picture, we use semi-analytical merger trees constructed using the 

■ extended Press-Schechter formalism, and evolve the subhalo populations using the 
^ ' average mass loss rates obtained from our simulations. The resulting subhalo mass 

• i-H , functions arc found to be in good agreement with the simulations. Our model can be 

' applied to semi-analytical methods of galaxy formation, to accurately follow the time 

H ' evolution of subhalo masses. 

^ ! 
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1 INTRODUCTION 

Understanding structure formation is a fundamental topic 
in modern cosmology. In the current ACDM concordance 
cosmology, the matter density of the Universe is dominated 
by cold dark matter (CDM), whose gravitational evolution 
gives rise to a population of virialized dark matter haloes 
spanning a wide range of masses. Numerical simulations of 
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structure formation in a CDM universe predict that these 
dark matter haloes contain a population of subhaloes, which 
are the remnants of halos accreted by the host, and which 
are eroded by the combined effects of gravitational heating 
and tidal stripping in the potential well of the main halo. 

Understanding the evolution of the subhalo mass func- 
tion, as function of cosmology, redshift, and host halo mass, 
is of paramount importance, with numerous applications. 
For one, subhaloes are believed to host satellite galaxies, 
which can thus be used as luminous tracers of the subhalo 
population. In particular, linking the observed abundances 
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of satellite galaxies to the expected abundance of subhaloes, 
provides useful insights into the physics of galax y forma- 
tion ( e.g.jMoore et al"] |_1999: B ullock et al. 2000: S omervilld 
l2002l : iKravtsov et alj 12006 : iVale fc Qstriken 120061 ). Stud- 
ies along these lines indicate that galaxy formation be- 
comes extremely inefficient in low mass haloes, and sug- 
gest that there may well be a large population of low mass 
subh a loes with no optical counterpart (e.g., Moore et al.l 



1 19991 : IStoehr et af] l2002l : IKravtsov et~al] l2004h . In princi- 
ple, though, these truly 'dark' subhaloes may potentially 
be detected via 7-ray emis sion due to dark matter annihila- 
tion in their central cores (IStoehr et al ]|2003l : lBerton3l2006l : 
?; IPieri et al.ll2007l : iDiemand et all l2007l ).~or via their im- 
pact on the flux-ra tio statistics of multiply-lensed quasar s 
(e.g., iMetcalf fc Madau 2001,: Da lai fc Kochanek I I200I ). 
Alternatively, these techniques may be used to constrain 
the abundance of subhaloes, which in turn has implica- 
tions for cosmological parameters and/or the nature of 
dark matter. The evolution of the subhalo mass function 
is also of imp ortance for the surviva l probability of disk 
galaxies (e.g., iToth fc Ostrikeil 1 19921 : iBenson.et al ] 12004 
[steward etall l2007h and even has implicati ons for direct 
detection experiments of dark matter (e.g., iGoerdt et al.l 
I2OO7I ) Finally, understanding the rate at which dark matter 
subhaloes loose mass has important implications for their 
dynamical fri ction times, and thus for the merger rates of 
galaxies (e.g.. IBenson.et all I2OO2I : IZentner fc Bullock! [ioO^ : 
iTavlor fc Babulll2004 ). 



Despite significant progress in the last years, there 
are still numerous issues that are insufficiently understood. 
What is the mass function of haloes accreted onto the main 
progenitor of a present day host halo? How do the orbits 
and masses of subhaloes evolve as they are subjected to 
dynamical friction, tidal forces and close encounters with 
other subhaloes? How does this depend on the properties of 
the host halo? In this paper we address these questions us- 
ing high-resolution numerical simulations. We trace back the 
evolution of self-bound substructures identified in present- 
day host haloes up to the point where they are first ac- 
creted by the main progenitor of the host halo. Using this 
method we are able to link the present-day population of 
subhaloes to the merging history of the host system. We will 
show that larger systems, forming at lower redshifts and so 
accreting their satellites more recently, co ntain at the end 
more subhaloes than small er hosts (see also lGao et al.|[2004l : 
Ivan den Bosch et alllioosh . 

In Section [5] we describe the simulations used. In Sec- 
tion [23] we present the algorithms employed to identify the 
haloes and to follow their merging history trees. In Section[3] 
we show how the unevolved subhalo mass function is con- 
structed from the merger tree of present-day halos and sug- 
gest an analytical fit. In Section [4] we explain how subhaloes 
are identified at the present time. In Section [S] we calculate 
the mass loss rate for subhaloes, and characterize its de- 
pendence on host halo mass and on redshift. In Section|6]we 
present Monte Carlo simulations that reproduce the subhalo 
mass function measured in the A''-body simulations. Finally, 
in Section [7] we summarize our results and draw some con- 
clusions. 



2 THE SIMULATIONS 

To study the evolution of the subhalo population of dark 
matter haloes we use two different types of simulations: a set 
of 48 massive haloes that have been extracted from a large 
cosmological simulation and resimulated at much higher res- 
olution, and a set of two large, cosmological simulations that 
probe a much larger dynamic range in host halo masses. 



2.1 Resimulations 

Our sample of 48 resimulated dark matter haloes is ex- 
tracted from ten high-resolution A''-body resimulations 
of galaxy clusters, containing 512'^ particles in a cube 
479Mpc//i on a side. All simulations assume a fiat ACDM 
model with flp = 0.3, h = 0.7, as = 0.9 and fib = 0.04 
(jYoshida. Sheth fc Diaferidl200l[ ). The masses of the haloes 
selected to be resimulated cover the range 5.1 X 10" - 2.3 X 
lO^'^ M0//1 at redshift z = 0. 

For the resimulations, we adopt a particle mass of 
1.3 X 10^ Mq /h and a gravitational softening length of e = 5 
kpc//i (Plummer equivalent). The initial conditions for the 
resimulation are generated with higher mass and force res- 
olution using the Z oomed Initial Condition technique (ZIC, 
iTormen etalll 19971 '): halo Lagrangian regions are populated 
with a larger number of less massive particles, and addi- 
tional small-scale power is added appropriately. The new 
initial conditions are evolved using the Tree-SPH code Gad- 
get2 (jSpringeli \200i) from redshift 2; = 60 to the present 
time, using dark matter particles only. We study these res- 
imulations using 88 output times equally spaced between 
2: = 10 and z — 0. Figure [1] shows an example of one of 
these resimulated cluster-sized haloes, embedded in its sur- 
rounding large-scale structure. The cluster is resolved with 
more than one million particles within its virial radius, and 
there are roughly 6 x 10^ hi gh resolution dark matter par- 
ticles inside 20 Mpc//i. See jPolag et al.1 |200i) for further 
details. 



2.2 Cosmological N-Body Simulations 

We will also make use of two cosmological A^'-body sim- 
ulations. The fir st is the so called "'GIF2"' simulation 
(jCao et all l2004h . a periodic cube of side 110 Mpc//i 
assuming the concordance ACDM model {0,a, ^bh^, h, 
o-g) = (0.7, 0.0196, 0.7, 0.9). The index of the initial power 
spectrum as be chosen n — 1, with the transfer func - 
tion produced by CMBFAST (jSeliak fc Zaldarriagal [l996l ). 
GIF2 contains 400"^ dark matter particles, each of mass 
1.73 X lO^Mo/Zi. We will use 50 output times logarith- 
mically spaced between z = 10 and z — 0, which suf- 
fices to construct accurate halo and subhalo m erging history 
trees using the method of lTormen et al.1 (j2004h . The numer- 
ical data of the GIF2 simulation are publicly available at 
http : //www. mpa-garching.mpg.de/Virgo 

Finally, as a consistency check, part of the analysis 
done on GIF2 was repeated on the lower resolution GIF 
box (ACDM each of mass 1.4 x 10^° Mq/K), which has the 
same cosmologica l para meters of the GIF2 simulation. See 
iKauffmann et aO (|l999[ ) for a detailed description of this 
simulation. 
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Figure 1. A cluster of the resimulated sample resolved with 2 X 10® particles within the virial radius and 6 X 10® high resolution particles 
within 40Mpc/h. 
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Table 1. The number of haloes in each logarithmic mass bin for the different simulations. For GIF & GIF2 we consider all haloes with 
a least more that 200 particles within their virial radius at redshift zero and whose main progenitor never has a virial mass exceeding 
the final value by more than ten percent. For the resimulated haloes we follow the merger tree and the satellites populations for all the 
haloes with more than 40000 particles at the present time. 

distance to the tenth closest neighbor. We assign to each par- 
ticle a local density pi.oM oc d~jo' ^'^^^ particles in density 
and take as centre of the first halo the position of the dens- 
est particle. We then grow a sphere of matter around this 
centre, and stop when the mean density within the sphere 



2.3 Halo-Finder & Merger History Tree 

We adopt the spherical overdensity criterion to identify 
haloes at each simulation output time (also called '"snap- 
shot"'). For each snapshot we estimate the local dark mat- 
ter density at the position of each particle by calculating the 
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falls below the virial value appropriate for the cosmological 
model at that redshift; f or the def i nition of virial density 
we adopted the model of I Eke et al.l (|l996l ): for example, at 
redshift = the virial density is pv = 324pi,, with pi, the 
mean background density of the universe. 

At this point we assign all particles within the sphere 
to the newly formed halo, and remove them from the global 
list. We take the centre of the next halo at the position of 
the densest particle among the remaining ones, and grow a 
second sphere. We continue in this manner until all particles 
are screened. We include in our catalogue only haloes with 
at least 10 particles within the virial radius; particles not 
ending up in any halo are considered as '"field"' or '"dust"' 
particles. 

We then build the merging history tree for all haloes in 
the simulation (or resimulation) using the halo catalogues at 
all snapshots, separated by redshift intervals dzi, as follows. 
Starting from each halo at z = 0, we define its progenitors 
at the previous output, z = dzi, as all haloes containing at 
least one particle that at z = will belong to that halo. The 
'main progenitor' ai z — dzi is defined as the progenitor 
that provided the largest mass contribution to the halo at 
z = 0. Next we repeat the same procedure, now starting 
at z — dzi and considering progenitors at z — dzi + dz2, 
and so on backward in time, always following the main pro- 
genitor halo. The resulting merger tree consists of a main 
trunk, which traces the main progenitor back in time, and 
of 'satellites', which are all the progenitors which, at any 
time, merge directly onto the main progenitor. 

In the following analysis we only consider haloes at red- 
shift z — whose main progenitor at any redshift has a virial 
mass M„(z) not exceeding the final value Mv{z = 0) = Afo 
by more than ten percent. In fact, any M^{z)/Mo signifi- 
cantly larger than one corresponds to an incomplete merger 
in which two haloes first merge but subsequently split again. 
This typically occurs when a (relatively) small halo pass 
through a larger one. Since these events complicate our anal- 
ysis, and their occurrence is only low, we decided to elimi- 
nate such merger histories. 

For all simulations we split the halo samples a,t z — 
in mass bins of width dlog(A^) = 0.5, with a minimum 
mass roughly corresponding to 200 particles within the virial 
radius for GIF and GIF2, and to 40000 particles for the 
resimulations. The actual mass bins for each run are listed 
in Table [U 



3 UNEVOLVED SUBHALO MASS FUNCTION 

Starting from each halo at z — 0, we trace its merger his- 
tory back in time and register all its satellites, i.e. all haloes 
directly accreted by the halo main progenitor at any out- 
put time. In order to remove subhaloes that at z = reside 
outside the host due to their elongated orbits (and so do 
not contribute to the subhalo population) we only consider 
satellites which donate at least 50% of their mass to the 
final halo. Let n{mv/Mo,z) be the number of satellites of 
virial mass mi,, accreted at redshift z by a host halo with 
mass Mo at redshift zero. Integrating this expression over 
the redshift interval zi ^ 2 ^ 22, we obtain the total num- 
ber of satellites of mass 7ti„ accreted by the main progenitor 



during that interval. 



N 



Mo 



(1) 



which we term the unevolved subhalo mass function. 

In Figure [5] we plot the unevolved subhalo mass func- 
tion for different redshift intervals, as measured in the GIF 
(left) and GIF2 (right) simulations. The data points refer to 
different mass bins of the parent haloes at redshift 2 = 0, 
as indicated. As stated above, we only considered satellites 
that contributed at least 50% of their mass to the final 
(2 = 0) host. Setting 21 = and 22 = 2max, which is the 
maximum redshift available in the simulations, we obtain 
the total, unevolved subhalo mass functions shown in the 
lower panels of Figure (5] Note that there is no indication for 
any significant dependence on Mo, indicating that the un- 
ev olved subhalo mass function is indeed universal, as found 
by Ivan den Bosch et al.l (|2005l ). After some experimenting, 
we find that the unevolved subhalo mass function is well 
fitted by 



dN 



d ln(m„/A/o^ 



= Nqx' 



aMo 1 



(2) 



with a = 0.8 and iVo = 0.21. This fitting function is indi- 
cated by solid black lines in the lower panels. 

In the upper and middle panels of Figure [2] we show 
the mass functions of the satellites accreted at redshifts 
larger and smaller than 2/, respectively. Here 2/ is the so- 
called assembly redshift, defined as the highest redshift at 
which the mass of the main progenitor Af„(2) exceeds half 
the final value, (i.e., M„(2) > A/o/2). Once again, the re- 
sults for different host halo masses are indistinguishable, 
and are extremely well fit by equation @ with a = 0.8. The 
normalization, A^o, however, needs to be adjusted. Naively 
one would expect that both mass functions should have a 
normalization No/2. However, because of the discreteness 
of the merger histories, the average mass at the forma- 
tion redshift. M(zf), is actually slightly larger than Mo/2. 
ISheth fc Tormeni \2004 ) have shown that, for the spherical 
collapse case and assuming a white noise power spectrum, 
the mass at formation has a distribution (eq.[4] of their pa- 
per) with a mean value usTOiMo = (0.586±0.005)Mo. Here, 
combining haloes from GIF and GIF2, we find a mean for- 
mation mass p,GiF+GiF2Mo = (0.572 ± 0.001)Mo. The fact 
that our results are somewhat lower owes to the fact that the 
dist ribution of p depends ( weakly) on the power spectrum 
(see lSheth fc Tormeni l2004l l. The normalizations of the un- 
evolved subhalo mass function accreted before Nq^i, and after 
A^o a the formation redshift are therefore: 



No.b = pNo = O.572A0 , 
No,a = (1 - p)No = O.428A0 . 



(3) 
(4) 



These are the normalizations that we have adopted for the 
fitting functions shown in the upper and middle panels of 
Figure H 



4 EVOLVED SUBHALO MASS FUNCTION 

The evolved subhalo mass function at any redshift 2 is built 
from all the satellite haloes accreted by the main progeni- 
tor at all redshifts larger than 2, where for each satellite we 
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Figure 2. Mass functions of accreted satellites (unevolved subhalo mass functions). In tlie panels the various data points and line types 
refer to different present-day host halo masses. In the figures the the bounds of the mass bins are expressed in unit of log(MQ//i). The 
solid lines represent the fitting function to the distributions: equation 1(2)1 (see the main text for more details). Note that we only consider 
subhaloes that at z = contribute at least 50% of their mass. This ensures that at z = their center of mass lies within the virial radius 
of the host. Top: Unevolved subhalo mass function accreted before the formation redshift 2/ of the host halo (defined as the earliest 
redshift when the mass of the main progenitor exceeds half the final mass). Center: Same as above, but only counting satellites accreted 
after Zf. Bottom: Same as above, but now counting satellites accreted at any redshift. 



compute at each redshift its self-bound mass msh{z)- Oper- 
ationally, we perform the following steps: 

• given a satellite halo, we identify its merging redshift, 
Zm, defined as the latest redshift when the satellite was still 
an isolated halo, just before it was accreted by the main 
progenitor; 

• we calculate the position of its cente r using the "mov- 
ing center method"' (jTormen et al.|[l997l ). i.e. by repeated 
calculation of its center of mass using smaller and smaller 
radii to identify the subhalo densest core; 

• we compute the subhalo tidal radius - as in 
iTormen et aTllj 19981 ): 

• we evaluate the binding energy of each subhalo particle 
by summing its potential energy (calculated using all par- 
ticles inside the tidal radius) and its kinetic energy (using 
its residual velocity with respect to the average value inside 
the tidal radius); 

• we remove all particles with positive binding energy, 
and iterate the previous steps until the self-bound subhalo 
mass converges. 



With these data in hand, we can follow the time evo- 
lution of the self bound mass of each subhalo, snapshot by 
snapshot, from the merging redshift Zm to the present time 
2 = 0. In the following Section we will use this informa- 
tion - gathered from the resimulated haloes - to estimate 
the subhalo mass-loss rates at all redshifts. 

In Figure [3] we show a schematic representation of the 
merging history tree of a halo. Time runs upward, with the 
final halo depicted at the top. Light gray circles indicate 
the main progenitor at each time, whose history defines the 
'main trunk' of the tree. Dark gray circles indicate satellite 
haloes, i.e. progenitor haloes accreted directly by the main 
trunk of the tree. The progenitors indicated by black circles 
are the 'leaves' of the merging-history tree, and reflect those 
progenitors whose mass is of the order of the resolution of 
the tree; i.e., for these haloes no progenitor can be identi- 
fied in the simulation at earlier outputs. Satellites are either 
'leaves', or they have progenitor haloes at earlier outputs, 
which in principle give rise to a population of sub-subhaloes, 
etc. In this paper we do not consider this substructure of 
subhaloes, though we intend to address their evolution in a 
forthcoming paper (Giocoli et al. in preparation). 
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Figure 3. Schematic representation of tlie merging-liistory-tree 
of an lialo. Solid light gray circles connected on the parent halo 
represent the main branch of the tree. Solid dark gray circles 
indicate satellites. Solid black circles indicate leaves progenitors. 
See the main text for explanation. 



As an example, Figure |4] shows the subhalo population 
of the most massive halo found at 2 = in the GIF2 cosmo- 
logical run. The left panel shows all particles inside the halo 
virial radius In the middle panel only those particles 
are shown that are bound to subhaloes located within J?„, 
while the right panel shows all other 'field' particles, which 
are bound to the main halo, but not to any subhalo. 

In Figure[^we plot the subhalo mass function for GIF2 
haloes at redshift z — 0, split according to the final halo 
mass. We considered all self-bound subhaloes with at least 
10 particles, and removed all subhaloes at halo-centric dis- 
tances r < 0.05i?„, where the subhaloes are difficult to de- 
tect. Note that the evolved subhalo mass function is not 
universal, but depends on the mass of the final host halo, 
with more massive haloes hosting more subhaloes of a given 
msb/Mo. 

The fact that the evolved subhalo mass function is not 
universal, but rather depends on host halo mass, Mo, was 
first shown by iGao et al] (|2004 ) . Using merger trees con- 
struc ted with the extended Press-Schechter (EPS) formalism 
fe.g.. lLacev fc Colelliggsl ). and adopti ng a simple model for 
the av erage mass loss rate of subhaloes. Ivan den Bosch et al.l 
l|2005l ) were able to reproduce these trends, which they ex- 
plained in terms of (i) the universality of the unevolved 
subhalo mass function, and (ii) the fact that more mas- 
sive haloes assemble later. The latter implies that smaller 
systems accrete their satellites at higher redshifts, when the 
haloes are denser, and the dynamical times are shorter. This, 
in turn, ensures that dynamical effects that promote mass 
loss are more efficient. Furthermore, satellites that are ac- 
creted earlier also are subjected to mass loss for a longer 
time. Both effects contribute to less massive host haloes hav- 
ing less substructure at 2; = 0. 

Below we will show that the subhalo populations in our 
resimulated haloes agree well with this picture, and that in- 
deed the average mass loss rates are higher at higher redshift. 
We will also find that the average mass loss rates are vir- 
tually independent of the mass ratio msb{z)/Mv{z) between 
the subhalo and its host. 



unevolved 




-0.5 



-3 -2 

» r / 1 

Figure 5. Subhaloes mass function of the self-bound particles of 
the haloes accreted by the main branch of the merger-history-tree 
of an halo, for GIF2 simulation. In the plot it has been considered 
all satellites with a distance from the center of the host halo less 
then the virial radius. We also plot the unevolved distribution: 
equation ((2]l. The different data points and line types used are 
the same of Figure |2] 



5 SUBHALO MASS LOSS RATES 

In this section we estimate the subhalo mass loss rate, mod- 
eling it as a function of (i) the instantaneous satellite to host 
mass ratio: msb{z)/Mv{z), (ii) the mass of the host halo at 
redshift zero: Mo, and (iii) the cosmic time through the red- 
shift z. For this purpose we will use the subhalo population 
identified in the resimulations, as haloes in this sample have 
better force, mass and especially time resolution (88 snap- 
shots between redshift ten and zero) than the cosmological 
GIF2 run. Since mass loss rates mostly depend on the local 
environment inside the host halo, the resimulated sample 
will provide correct rates even if the haloes themselves do 
not necessarily represent a fair sample for the given cosmo- 
logical model. 

In Figure |H] we show the unevolved subhalo mass func- 
tion for satellites identified in the merger trees of the set 
of resimulated haloes; host haloes are split in three mass 
bins, according to Table [T] As for the GIF2 simulation, the 
unevolved subhalo mass function obtained from the resimu- 
lations is well fit by eq. [S] 

After a satellite enters the virial radius of the host, var- 
ious dynamical effects, including dynamical friction, tidal 
stripping, and close encounters with other subhaloes, cause 
the subhaloes to loose mass, a nd may eventuall y result in 
their complete disruption (e.g.. IChoi et al.|[2007[ ). The (av- 
erage) mass loss rate of dark matter subhaloes is the di- 
rect link between the unevolved and evolved subhalo mass 
functions, and also is a fundamental ingredient for semi- 
analytical models of galaxy formation, as it sets the rate 
at which satellite galaxies merge with the central galaxy in 
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Figure 4. Subhalo population. Left: all particles composing tiie most massive halo found at 2 = in the GIF2 simulation; the virial 
mass for this halo is A'ly = 1.8 X 10^^ Mq/Zi, resolved by more than one million particles. Center: particles bound to subhaloes at redshift 
z = 0. Right: particles bound to the main halo but not to subhaloes. 



I ' ' ' ' I 

__13,7 - 14 




Figure 6. Unevolved subhalo mass function for the resimulated 
haloes. We notice that the function is independent on mass and 
well described by the same function fitting the GIF2 data (Fig- 
urc[2}. Haloes are split in three mass bins. In the figure the bounds 
of the bin are expressed in unit of log(MQ/h). 



a halo, it determines the evolution of the mass-to-light ra- 
tios of satellite galaxies, and it regulates the importance of 
stellar streams in the haloes of central galaxies. 

In this section we measure the mass loss rate experi- 
enced by each satellite. In addition, using statistical averag- 
ing, we determine the average mass loss rate of satellites as 
a function of the parameters listed at the beginning of this 
section. We define the average mass loss rate between two 



d / msb\, 



successive snapshots at redshift, zi and Z2, as 
msb(z2) _ msbjzi) 

j-^ ; , Zi <Z<Z2. 5 

t{Z2) - t[Zi) 

where the values of mstiz) and Mv{z) axe obtained by linear 
interpolation of the values at Zi and Z2. In Figure [7] we 
plot the subhalo mass loss rate as a function of the ratio 
msb{z)/Mv{z); Each panel refers to a different bin for the 
mass My (z) of the host halo. 

The green points and band in each panel indicate the 
median and quartiles of the distribution. The thick magenta 
line represents a least squares fit to the median values in each 
panel; the fit is limited to the region where the median ex- 
hibits a linear behavior; we excluded by hand median values 
for rUsb/Mv close to one, which correspond to major merg- 
ers and cannot be described by a simple mass loss model. 
The thin dashed black line, identical in all panels, shows 
the global least square fit obtained using the data from all 
panels. 

The data show a clear linear relation between rUsb/Mv 
and its time derivative, so we can write our model as; 



log 



d{msb/Mv) 



dt 



■ alog{msb/My) + b. 



(6) 



Exponentiating this relation, and expanding the derivative 
on the LHS, we obtain; 



msb _ My nisb 

My My My 



= 10° 



My 



(7) 



Due to the large number of snapshots in the resimulations, 
the time separation between two subsequent snapshots is 
always short; dt ~ 0.1 Gyr. This is small enough to assume 
a constant mass for the host halo; My = 0. By doing so, we 
obtain an expression for the specific mass loss rate 

(8) 



msb 
rrisb 



msb 

My 



where the free parameters T{z,My) — 10"* and ^(0,M„) = 
a — 1 might in principle depend both on cosmic time (or, 
equivalently, redshift z) and on the virial mass A{y{z) of the 
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Figure 7. Subhalo mass loss rate. Each panel refers to a different bin in host halo virial mass at the redshift when the mass loss rate is 
computed. The filled circles represent the median of points and the hatched region the quartiles. The thick solid line is the least square 
fit to the median distribution for each panel. The thin dashed line is the average least squares for the different host halo masses. 



host halo at that time. The negative sign arises from the 
mass loss of the satellites. Not e that this specific rnass los s 
rate is identical to that used bv lvan den Bosch et al.l (|2005h . 

Figure [8] shows how the time scale of the mass loss rate, 
r = 10~\ and ( = a — 1, as measured from the data shown 
in Figure [T] depend on the virial mass, , of the instanta- 
neous host halo. Error bars reflect the usual uncertainty on 
the coefficients obtained from the least square fitting. The 
slope is found to be independent of the mass of the host halo, 
with a best fit value of a = 1.07 ± 0.03 (( = 0.07 ± 0.03). 
This implies that the specific mass loss rate is almost inde- 
pendent of the instantaneous mass ratio rrist/Mv. On the 
other hand, the zero point, b, is found to be larger for less 
massive haloes. 

In order to show the typical spread of points in each 
panel of Fig. [7] around each median, in the bottom panel of 
Figure |5] we show the average (over the six panels of Figure 
O of the differences between each quartile and the median 
itself. We see that on average fifty percent of the points lay 
roughly within a distance log y = ±0.3 from the median; 
that is, typical mass losses deviate from their median value 
by less than a factor of two. 

In Figure[9]we plot the subhalo mass loss rate versus the 
ratio msb{z)/Mv{z), now binned according to the redshift 



at which the mass loss rate is calculated. Medians, quartiles 
and lines are as in Figure [T] The time scale r and C, for 
the six panels are shown in Figure 1101 plotted versus the 
mean redshift of each of the six bins; in the bottom panel 
the average quartile distribution for each fit (as explained 
above) is shown. The red solid curve superimposed to the 
trend in zero point is the equation 



t{z) = To 





-1/2 


\H{z)] 


-1 


[ Ao J 




[ Ho \ 





(9) 



with H{z) the Hubble constant at redshift z, and with to = 
2.0 Gyr. 

This equation was proposed by Ivan den Bosch et al.l 

(|2Q05h and describes the redshift dependence of mass loss 
rates obtained under the assumption that r is proportional 
to the dynamical time tdyn oc p„ ^^^(z), taking into account 
that, according to the spherical collapse model, the aver- 
age density within the virial radius, pv is independent of 
halo mass at fixed redshift. This means that we can write 
t{Mv,z) = t{z). The red line in Figure (TU] shows that in- 
deed this provides a good description of the measured mass 
loss rates. 

Note, though, that Figure [5] suggests that the average 
mass loss rates also depend on host halo mass. In order to 
reconcile this with the claim that the zero-point is indepen- 
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Figure 9. Subhalo mass loss rate. Each panel refers to a different bin in the redshift at which the mass loss rate is computed. Symbols 
and lines are as in Figure [7] 



dent of M„, recall that, on average, more massive haloes 
assemble (and thus accrete their satellites) earlier than less 
massive haloes. Therefore, the different panels of Figure [7] 
actually refer to different average redshifts, with larger 
corresponding to a lower average redshift. Consequently, the 
'apparent' mass dependence evident in the upper panel of 
Figure |S] is merely a reflection of the redshift dependence 
described by equation (|9]). To demonstrate this we now split 
the data points of each panel of Figure[9]in different subsets, 
according to the mass of the host halo. Figure [TT] shows the 
average slopes and zero points obtained for these subsets 
using least-squares fitting. This clearly shows that the char- 
acteristic time scale for mass loss (given by the zero point) 
is independent of the host mass Mv{z) at fixed redshift, in 
accord with equation ([9|. 

Thus, to good approximation, the average mass loss rate 
of dark matter subhaloes depends only on the density of the 
host halo, and thus on redshift (or cosmic time), but not on 
the mass of the host halo. Furthermore, since the best-fit 
value of C, is close to zero, to good approximation subhalo 
masses decay exponentialljQ according to 



'nisbi't) = iTiv exp 



(10) 



^ this follows from a simple integration of equation ^ with C = 



where m„ is the mass of the satellite at the time of accretion, 
tm, and t[z) is given by equation @ with tq = 2.0 Gyr. 



6 MONTE CARLO SIMULATION 

In this section we compare our results to those of 
Ivan den Bosch et all (j2005h . and we use their Monte Carlo 
method to check the self-consistency of the results presented 
above, i.e., we check whether the (universal) unevolved sub- 
halo mass function, combined with the satellite accretion 
times and the average mass loss rates, can reproduce the 
evolved subhalo mass functions presen ted in Section [H 

The Monte-Carlo method of Ivan den Bosch et al.l 

()2005l ) starts by constru cting EPS me r ger t rees using 
the method described in [ van den BoschI (|2002l ) (see also 
ISommerville fc Koiatdliggg T These merger trees are then 
used to register the accretion times and masses of satel- 
lites m erging onto the main progen itor. Starting from these 
inputs, Ivan den Bosch et al.l (|2005l ') then proceeded as fol- 
lows. In between two time-steps, they evolve the masses of 
the subhaloes using equations [H and |9] The two free pa- 
rameters. To and C, were tuned to reproduce the subhalo 
mass function of massive, cluster sized halo e s obtained from 
numer ical simulations by ICao et af] (|2004h . IPe Lucia et al.l 
120041 ') and lTormen et al.l (|2004l 'l. This resuhed in ro = 0.13 
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Figure 8. Dependence of the fit parameters of the Figure [7] on 
the host halo virial mass. The top panel shows the time scale of 
the mass loss rate r = lO"*". The average and the least square fit 
of the data points have been computed on the plane (b, M^). In 
the central panel we show the dependence of parameter C, = a — 1 
on My. In the bottom panel we show the spread of the first and 
third quartiles around the median, averaged over the six panels 
of Figure [7] (see the main text for a detailed explanation) . 
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Figure 10. Time scale of the mass loss rate and C, in term of 
the rcdshift at which the subhaloes are experiencing mass loss 
(Figure |9] The average and the least squares fit of the top panel 
were computed on the plane (b, z). The bottom panel shows the 
average first and third quartile for the median distribution in each 
panel of the Figure O constructed as previously described in the 
main text. 
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Figure 11. Time scale of the mass loss rate and f versus host 
mass, for six fixed redshift bins - represented by different data 
points. The horizontal lines, with various line type, show the av- 
erage 6 = — log(T) and C, for each redshift bin. 



Gyr and C, = 0.36, which diflers substantially from the re- 
sults obtained here: ro = 2.0 Gyr and = 0.06. The reason 
for this discrepancy owes to the use of EPS merger trees, as 
opposed to merger trees extracted from numerical simula- 
tio ns. In fact, the unevo l ved su bhalo mass function obtained 
bv lvan den Bosch et al.l (|2005l l differs significantly from that 
shown in Figures [5] and [S] in that it is significantly higher, 
and with a different slope at the low mass enqj. Conse- 
quently, in order to reproduce the su bhalo mass functions 
obtain ed from numerical simulations, Ivan den Bosch et al.l 
llooi) had to adopt higher mass loss rates (i.e., a lower 
value for ro, and a different mass dependence (i.e., a differ- 
ent value for C^. 

The fact that EPS merger trees predict an unevolved 
subhalo mass function that differs significantly from that 
obtained in numerical simulations, should not come en- 
tirely as a surprise. After all, the construction of EPS 
merger trees relies on the spherical collapse model (see 
iLacev fc Colell 19931 : [Sommerville fc Kolattlll999l ). However, 
in reality, the collapse of dark matter haloes is influenced 
by the surrounding tidal force field, ma king the collapse 
ellipsoidal, rather than spherica l (e.g., Sheth fc TormenI 
1 19991 : ISheth. Mo fc TormenI [2OO1I: ISheth fc TormenI I2OO2I ). 



As shown bv ISheth fc TormenI lj2002h ~ the conditional and 
unconditional mass functions are different under ellipsoidal 
collapse conditions than under spherical collapse condi- 
tions, which has important consequences for the accuracy 
of the EPS merger trees. For instance, the halo formation 
times predicted by EPS are systematically offset from those 



^ It is noteworthy, though, that the EPS formalism predicts that 
the unevolved subhalo mass function is universal, i.e., indepen- 
dent of the host mass, in good agreement with the simulation 
results presented here. 
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Figure 12. The dotted histogram sliow tlie mass accreted by tiie 
main branch in the Monte Carlo merger tree with the overplotted 
equation l(2}. The sohd Unes represent the subhaloes mass func- 
tion obtained evolving the mass accreted by the main progenitors 
of different present day Mo-halo. 



obtained from numerical simulations (iLacev fc Cole 1 1 19931: 
Somerville et al.ll2000l: Ivan den Boschll2002l : IWechsler et all 



20021 : iGiocoli et all \200'K while the aver age mass of the 



main progenitor is typically overestimated (jSomerville et all 
l200d ). 

To perform the self-consistency check mentioned 
ab ove, we therefore us e the same Monte-Carlo method 
as Ivan den Bosch et al.l lj2005h . but we randomly remove 
satellite-branches from the merger tree with a probability 



reject — 



nEPs{mv/Mo) 



(11) 



where nsim and tieps are the unevolved subhalo mass func- 
tions obtained from the simulations and from the EPS 
merger trees, respectively. This ensures that the Monte 
Carlo method uses an effective, unevolved subhalo mass 
function t hat is identical to that of eq uation 

As in Ivan den Bosch et al.1 |200i) we evolve the masses 
of the subhaloes using equations [8] and |9] with to = 2.0 Gyr 
and = 0.06, which are the best-fit values obtained in sec- 
tionO The resulting evolved subhalo mass functions, for five 
different masses of the present-day host halo, are shown in 
Figure 1121 together with the unevolved subhalo mass func- 
tion obtained using the rejection scheme outlined above (and 
which is independent of the host halo mass). Each evolved 
subhalo mass function is th e average obtained from 2 000 
merger tree realizations (see Ivan den Bosch et al]|2005l . for 
details). A comparison with the evolved subhalo mass func- 
tions obtained from our numerical simulations, and shown 
in Figure [T^ shows good agreement. This indicates that the 
evolved subhalo mass functions are self-consistent with the 
(universal) unevolved subhalo mass function and the simple 
form for the average mass loss rate obtained in this paper. 



7 SUMMARY AND CONCLUSIONS 

In this paper we have studied the mass loss rate of dark 
matter subhaloes using a set of high resolution A'^-body sim- 
ulation of structure formation. Haloes were followed back- 
ward in time along the main branch of their merging history 
tree. At each snapshot the satellites accreted by the main 
branch were identified. We showed that the mass function 
of accreted satehites (unevolved subhalo mass function) is 
universal, that is, it does not depend on the present day host 
halo mass Mo, and we presented a fitting function for this 
distribution. 

We then followed each accreted satellite forward in time, 
snapshot by snapshot, computing its self-bound mass and its 
mass loss rate. We fo und that the express i on for the mass 
loss rate proposed bv Ivan den Bosch et al.1 (120051 ) IS consis- 
tent with A^-body simulations, and excellent agreement is 
obtained if the value with tq = 2.0 Gyr is taken. In addi- 
tion, we find that the average mass loss rate is virtually 
independent of the instantaneous mass ratio rUsb/My be- 
tween the subhalo and its host halo. This differs substan- 
tially from the best-fit mass loss rate parameters obtained by 
Ivan deii Bosch et al.l (120051') using E PS merger trees. In par- 
ticular, Ivan den Bosch et al.1 (|2005l ') obtained tq = 0.13 Gyr, 
and a significant dependence on msb/Mv The reason for this 
discrepancy is that the unevolved subhalo mass function of 
EPS merger trees is too high, so that a higher mass loss rate 
was inferred to be consistent with the evolved subhalo mass 
functions in numerical simulations. 

With an unevolved subhalo mass function that is uni- 
versal, and an average mass loss rate that is virtually inde- 
pendent of nisb/My, it becomes straightforward to under- 
stand why less massive haloes have evolved subhalo mass 
functions with a lower normalization. This simply owes to 
the fact that less massive haloes assemble earlier, which im- 
plies that they accrete their satellites earlier. At earlier times 
the mass loss rate is higher, because the dynamical times of 
dark matter haloes are shorter. In addition, a subhalo that is 
accreted earlier is subjected to mass loss for a longer period. 
Both these effects contribute to the fact that less massive 
haloes have less substructure. 

The present description does not consider the possible 
presence of subhaloes within subhaloes, accreted along the 
tree of each present day subhalo. In a follow-up paper (Gio- 
coli et al. 2008, in preparation) we will investigate this issue 
in detail, comparing the populations of subhaloes found us- 
ing different techniques. 
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